Driven Diffusion in Periodic Potentials with Stochastic Path Integral Hyperdynamics 
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We consider the driven diffusion of Brownian particles in ID periodic potentials using the recently 
proposed Stochastic Path Integral Hyperdynamics (SPHD) scheme [L.Y. Chen and L.J.M. Horing, 
J. Chem. Phys. 126, 224103 (2007)]. First, we consider the case where a single Brownian particle 
is moving in a spatially periodic potential and subjected to an external ac driving force. We confirm 
that there is no stochastic resonance in this system and find that at higher frequencies the diffusion 
coefficient D is strongly suppressed. The second case is that of a dimer moving in a periodic potential 
with a static bias. For this case, there's a strong suppression of D when the dimer bond length is 
an integer multiple of the lattice constant of the potential. For both cases, we demonstrate how the 
SPHD allows us to extract the dynamical information exactly at different bias levels from a single 
simulation run, by calculating the corresponding statistical re-weighting factors. 

PACS numbers: 05.10.-a, 05.40.Jc, 05.10.Gg, 87.15.Vv 



INTRODUCTION 

The study of particles performing Brownian motion 
in a periodic potential constitutes a hallmark example of 
stochastic particle dynamics with important applications 
in various branches of science and technology. Perhaps 
the most common application of periodic Brownian mo- 
tion is the diffusive dynamics of atoms and molecules on 
crystal surfaces Surface diffusion is among the most 
important mechanisms that controls processes such as 
island nucleation and subsequent surface growth. It has 
been shown that by controlling the mobility of particles 
on the surface by external means, such as an ac or dc elec- 
tric field, allows morphological control over the growing 
surfaces [lj . It is thus of great interest to model periodic 
Brownian motion with static and time-dependent exter- 
nal fields. 

To this end, there have been several studies reporting 
the diffusion of a single Brownian particle in a periodic 
potential with external ac bias applied [111 HQ. The 
case of a dimer consistingof two connected particles has 
also been considered @, 0, @|, in which case there's an 
additional length scale in the problem, namely the dimer 
bond length, as shown schematically in Fig. 1. Most 
of the studies reporting the behavior of Brownian parti- 
cles discuss the influence of an oscillating bias on trans- 
port coefficients. The central issue here is existence of 
a stochastic resonance (SR), which can greatly enhance 
the diffusion coefficient D in 2D Q. However, it has been 
shown in the case of ID periodic potentials that although 
the local jump rate of particles can be enhanced, there is 
no true SR in the hydrodynamic limit 0, [B[ . An interest- 
ing limit of the periodic Brownian motion is where the 
energy barrier Vq is much larger than the thermal energy, 
i.e. W » 1 0, where = l/k B T and k B is the Boltz- 




FIG. 1: Driven diffusion in a ID periodic potential. The rele- 
vant parameters are the barrier height Vo, the lattice constant 
A and the dimer bond length ao. The details of an external 
bias are discussed in the text. 



mann constant and T the temperature (cf. Fig. 1). Since 
Brownian motion is activated by thermal fluctuations, 
the diffusion rate is proportional to exp(— /3Vb) which be- 
comes very small at low temperatures. To overcome this 
rare event problem in Molecular Dynamics (MD) sim- 
ulations, Voter [9( has proposed the so-called Hyperdy- 
namics (HD) scheme, which involves accelerating the dy- 
namics by adding proper bias potential, which effectively 
lowers the barrier height. The dynamics is then cor- 
rected based on the approximate Transition State Theory 
(TST). There exist various approaches to the choice of 
the bi as pote ntial, and some examples can be found in 
Refs. 



bi as p otential, ai 
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However, recently a new scheme has been proposed 
that is based on the mapping of the stochastic Langevin 
equation to a path integral form [Io| . Unlike the standard 
HD scheme, this so-called Stochastic Path Integral Hy- 
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perdynamics (SPHD) method allows an exact correction 
of the dynamics by resampling the simulated paths. In 
other words, this method is not restricted to the TST ap- 
proximation. Further, it is not restricted to static energy 
barriers; both entropic barriers and even time-dependent 
bias can be employed. This allows an efficient way to 
overcome the large barrier problem, as demonstrated in 
Ref. [1(3]. However, the SPHD method is not limited to 
the case of high barriers. Since in principle any external 
bias force can be used, it should be possible to obtain 
results for many different bias values from running LD 
simulations with a single value of the bias force, or even 
without such a force if need be. To demonstrate this, 
in this work we have undertaken to employ the SPHD 
method to study periodic Brownian motion for two in- 
teresting cases. In the first case, we consider the diffusion 
of a single Brownian particle in a ID periodic potential 
with ac forcing. The second case is that of a dimer diffus- 
ing in a ID periodic potential with a static bias. For both 
cases, we show how the SPHD method can be efficiently 
employed to obtain the transport coefficients for a range 
of different external forcing terms from single simulation 
runs. 

STOCHASTIC PATH INTEGRAL 
HYPERDYNAMICS 

Brownian motion of a single particle can be represented 
by the Langevin equation 

mr(t) + mjr(t) - F = £(t), (1) 

where r(t) denotes the position of the particle (of unit 
mass m) at time t, moving under the influence of ex- 
ternal force F and the random noise term £(t) satisfies 
(£(*)) = and (£(*)£(*')) = 2k B TmjS(t - t'), where 7 
denotes the friction coefficient. In general, to accelerate 
the dynamics one must add a bias force Fb(r,t) to the 
Langevin equation to get 

mr(t)+mjr(t)-F-F b (r,t)=£_{t). (2) 

With an appropriately chosen bias potential, the dynam- 
ics of the system evolves much faster than in the original 
Langevin equation. The dynamics given by numerically 
solving Eq. (2) is fictitious, of course. However, the use 
of path integral formalism allows an exact compensation 
of the effect of adding the bias force by defining an effec- 
tive action functional I^(t) [13] 

i 

This expression constitutes time integration for a given 
realization of the random noise force £ along a given tra- 
jectory r(t). To recover true dynamics in the absence of 



Fb(r,t), one has to estimate the SPHD statistical weight 
factor exp(— /3if) and simply use it to re- weight every 
sampled trajectory. 

The fundamental quantity associated with Brownian 
dynamics is the single-particle (tracer) diffusion coeffi- 
cient which can be defined through the mean square 
displacement (MSD) of the tracer particle as 

D= lim l([r(t)-r(0)] 2 ). (4) 

t^oo Zt 

When studying particle diffusion the mean square dis- 
placement (MSD) at zero bias (true dynamics) can be 
obtained by running SPHD with a bias, calculating 
along every trajectory and re-weighting as 

([r(0) - r(t)]l) = ([r(0) - r(t)] 2 HD e~^) , (5) 

where the subscripts and HD correspond to the quan- 
tities with zero bias (true dynamics) and finite bias, re- 
spectively. The SPHD has an additional powerful feature 
which can be seen from the biased Langevin equation and 
the expression for the re- weighting factor 1^ (i) . It is pos- 
sible to obtain the true dynamics for any value of the bias 
force < fi, < Ff, by estimating I^(t, fb) corresponding 
to some fb- With this, the true dynamics at multiple bias 
levels (fb) can be extracted, simultaneously, from a sin- 
gle simulation run of biased dynamics according to Eq. 
(2) and estimating the weight factors, corresponding to 
multiple values of fb along all the trajectories. 

RESULTS AND DISCUSSION 

Brownian Particle in Periodic Potential with 
Time- Varying Bias 

The first case where we consider the application of the 
SPHD method is that of a Brownian particle in a one- 
dimensional spatially periodic potential with an external, 
time-dependent ac driving force 0, H, fl Hi • For such a 
system, the equation of motion is given by 

mx(t) + mix{t) — F = £(t) + Asin(2inot), (6) 

where the second term on the right hand side indicates 
an ac driving force with amplitude A and frequency 
lu. The diffusion of a Brownian particle can be stud- 
ied with respect to various values of these two param- 
eters. F is force due to spatially periodic potential 
— (Vb/2)[1 — cos(27nc/A)]. Here, we have employed the 
SPHD method by numerically solving Eq. (6) with A = 
using the standard velocity Verlet scheme |15j. The dif- 
fusion coefficients for different values of A and lu can be 
obtained by choosing fb(A, w) = A sin(27ra;i) and then es- 
timating the correction factor I^(t, f b ) and the reweighing 
factor exp(— 01$) for every bias force. 
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FIG. 2: Diffusion coefficients for the ac-driven Brownian par- 
ticle in a periodic potential as a function of the driving fre- 
quency u for various values of the amplitudes A. The hori- 
zontal dashed line shows the reference value of D for A = 0, 
and the vertical dashed line indicates the position of the res- 
onance frequency c3 r (see text for details). The solid lines are 
just guide to an eye and contain no other significance. 



We set scales for length = A, energy = kgT and mass = 
m. Then time unit is defined as t = A \J (m/fcsT) and all 
other relevant quantities are expressed as dimensionless 
and indicated by a tilde over the respective symbols. 

The parameters we have used in the present work are 
Vq = 2, and 7 = 2. The time step we have chosen is 
At = 0.0005 [16j. From the unbiased runs (in this case 
with A = 0) we have computed the MSD and the corre- 
sponding diffusion coefficients as given in Eq. (4) for a 
range of values of A (A = 0.1,0.5,1.0,1.5) and the fre- 
quency, as summarized in Fig. 2. All the data have been 
averaged over 10 6 trajectories. Our results demonstrate 
that for a given amplitude A of the external ac driving 
force, D decreases monotonically with the frequency ui. 
For higher values of u) studied, it appears that large val- 
ues of amplitude A lead to a lower D than small values 
of A. It seems that the ac bias for higher values of both 
A and u) acts detrimental to activation in the diffusive 
motion. This feature can be effectively used to localize 
the motion of the particle. 

It is of interest to compare the SPHD with direct sim- 
ulations for the Langevin equation. In this case, a larger 
time step of At = 0.005 was sufficient, and averages were 
taken over 10 6 trajectories. In Table I we show results 
for the case of Q = 0.2 as a function of increasing ampli- 
tude A. We find excellent agreement between these data. 
For the highest value of A, studied here, it appears that 
SPHD gives a slightly smaller estimate of D than the 
LE, although the data still agree within the error bars. 
Our preliminary results indicate that with increasing bias 
path sampling should be increased, too, which is caused 
by the exponential decrease in the reweighing factor with 
increasing A. 



TABLE I: Values of the diffusion coefficient D from SPHD 
and direct solution of Eq. (6) with ui = 0.2. 



A 


D 


D 




(from Eq. (6)) 


(SPHD) 


0.1 


0.155 ±0.001 


0.153 ±0.002 


0.5 


0.163 ±0.002 


0.159 ±0.004 


1.0 


0.185 ±0.003 


0.178 ±0.005 


1.5 


0.225 ± 0.003 


0.215 ±0.007 



An interesting issue in Brownian motion under time- 
periodic forcing concerns the existence of stochastic res- 
onance, which leads to a significant enhancement of the 
relevant transition rates In the case of a double- well 
potential, SR is expected to occur in the vicinity of the 
matching condition uj r — irr e , where r e is the (thermal) 
escape rate @. In the case of an extended periodic po- 
tential there's enhanced escape for local diffusion jumps 
over the barrier Vq 0. However, it has been shown in 
Refs. [3, H, @| that this enhancement exactly cancels out 
in the hydrodynamic limit for a ID periodic potential 
such as used in the_present study. One can estimate the 
frequency (u> r ) as 



e_p 
I 
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(7) 



The calculation results in uj r = 1.336. As mentioned 
above, there is no peak in the value of D around ui r which 
confirms the absence of stochastic resonance. 



Brownian Dimer in Tilted Periodic Potentials 

The second case that we consider here is that of a 
Brownian dimer diffusing along a tilted ID periodic po- 
tential (cf. Fig. 1). The equations of motion for the 
beads are given by 

mXi(t)+myxi(t)-F = &(t)-V x< [—{l-cos(2nXi/ \))-b t x 

(8) 

where i = 1,2 and the variables cc,(t) denote the posi- 
tions of the two beads. The second term on right hand 
side denotes the force associated with external static bias 
with spatial periodicity A (see Fig. 1) and tilt bf . On the 
left hand side, F represents the force due to the inter- 
action between the constituent monomers of the dimer 
which is given by a combination of the Lennard- Jones 
(LJ) potential and the FENE (Finitely Extensible Non- 
linear Elastic) potential. The LJ potential is given by 



U LJ (r) = 4e 
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r J \r 
0, otherwise. 



e, for r < 2 1 / 6 cr; 



(9) 



4 



0.05 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 



FIG. 3: The dimer diffusion coefficient D vs. the wavelength 
do /A for various tilts. Error bars are of the order of the size 
of the symbols, or smaller. The solid lines are just guide to 
an eye and contain no other significance. 



Here, e and a define the energy and length scales, respec- 
tively, and r is the separation between the monomers. 
The FENE potential is given by 



Ufene{t) 



1 



kR 2 Hl-r 2 /R 2 ), 



(10) 



where k is the effective spring constant and i?o sets the 
maximum allowed separation between the monomers. 

The SPHD method was employed by numerically in- 
tegrating Eq. (8) corresponding to a spatially peri- 
odic surface with given A and zero tilt. The dynamics 
along tilted periodic potentials was recovered by setting 
fb(bt) = — V Xi [— btxi\. For the resampling, the corre- 
sponding weight factors were evaluated as in the previous 
section. 

Again, we set scales for length = a, energy = s and 
mass = m. Then time unit is defined as t — a^J (jn/e) 
and all other relevant quantities are expressed as dimen- 
sionless and indicated by a tilde over the respective sym- 
bols. 

The parameters were set tobeV^=T = 0.1,7 = l 
and the time step At — 0.0005, with averages taken over 
10 6 trajectories. The FENE parameters are k — 1 and 
i?o = 2. The periodicity of the ID potential was changed 
between A = 1/2 and 2, and the tilt parameters were 
bt = 0.0,0.05,0.1. The effective diffusion coefficient was 
estimated by using the relation 



D = lim 



(x 2 



(*)) 



2t 



(11) 



where x CTn is the center of mass of the dimer, in order 
to subtract the drift term caused by b t > 0. In Fig. 
3 we summarize our results for D as a function of the 
wavelength of the periodic potential for three tilt val- 
ues. For the case of a dimer, there are now two relevant 
length scales in the system (see Fig. 1): the zero tem- 
perature equilibrium bond length of the dimer oq ~ 1.10 



and the wavelength A of the underlying periodic poten- 
tial. Thus, for the dimer motion there's a matching of the 
two lengths when the ratio ao/A = n is an integer. On 
the other hand, when this ratio is a half-integer, there's 
strong competition between the dimer bond and poten- 
tial energy. In the latter case, it should be easier for the 
dimer to escape as the effective diffusion barrier is lower. 
Indeed, as seen in Fig. 3 we find that for every value of 
the tilt there is a strong decrease in D in the vicinity of 
the first matching condition ao/A = 1 [17]. On the other 
hand, near ao/A = 3/2 there is a local maximum in D, 
as expected, followed by another pronounced minimum 
near ao/A = 2. In our model, the interaction between 
the monomers is highly anharmonic. This results in an 
equilibrium dimer separation which is temperature de- 
pendent and the first minimum in D in Fig. 3 is shifted 
from the zero temperature value ao/A = 1 to ao/A = 1.2. 
A finite tilt increases the overall magnitude of D while 
its non-monotonic behavior as a function of ao/A prevails 
almost unchanged . Our results are consistent with the 
Langevin dynamics studies of Bammert et al. [g], who 
considered dimer diffusion in a 2D periodic square poten- 
tial (with a hydrodynamic interaction term included) and 
found that there's a local maximum in D for ao/A = 3/2. 
Heinsalu et al. Q have reported dimer diffusion on a ID 
washboard-like potential, and they also find a flat mini- 
mum at ao/A = 1 for small tilts. 

Another interesting feature in our results is that the 
minima in D deepen with increasing n. This can be un- 
derstood as follows. When the matching condition is met 
(i.e. n is an integer), both beads in the dimer can sit ex- 
actly at the minima of the potential separated by n — 1 
minima. The elementary diffusion move of the dimer (lo- 
cal jump rate) consists of both of the beads crossing the 
saddle points (potential maxima) synchronously (assum- 
ing that bond fluctuations can be neglected). For every 
n the effective diffusion barrier is exactly the same in- 
dependent of n. However, the jump length of the dimer 
depends on n (in units of a~o) such that the distance the 
dimer moves is given by do/n. Thus, if we use the Dy- 
namical Mean Field theory to approximate the diffusion 
coefficient 18], the prefactor of D is proportional to the 
jump length squared, which gives 



D cx \. 



(12) 



We have analysed the data upto around ao/A = 4 and 
checked that minima of D near integer values of ao/A 
indeed decreases with increasing n, but slower than pre- 
dicted by Eq. (12) for the present set of parameters. 
Indeed, Eq. (12) strictly holds only in the limit of a rigid 
dimer bond and Vq 3> 1, which is not in the range of the 
parameters used here. 

In analogy to the single particle case, we have com- 
pared our results from the SPHD scheme to data ob- 
tained from directly integrating the Langevin equation 
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TABLE II: Comparison between results for the diffusion co- 
efficient D obtained from SPHD and from Eq. (8) with tilt 
b t = 0.10. 



a /A 


D 

(Langevin eq.) 


D 

(SPHD) 


1/2 


0.0482 ± 0.0004 


0.0470 ± 0.0005 


2/3 


0.0454 ± 0.0005 


0.0450 ± 0.0004 


1 


0.0368 ± 0.0003 


0.0370 ± 0.0004 


3/2 


0.0393 ± 0.0004 


0.039f ± 0.0006 


2 


0.0320 ± 0.0004 


0.0323 ± 0.0006 



with the static bias for the largest value of the tilt 
b t = 0.10, for various values of A. For the straight LD 
simulations At = 0.005 and averages were taken over 
10 6 trajectories. There is again good agreement between 
these two sets of data. 

An intriguing further extension of the SPHD method is 
extrapolation with more than one bias force parameter. 
In the present case, we did some test runs by numer- 
ically integrating Eq. (8) corresponding to a flat sur- 
face and zero static bias (Vo = b t = 0). The dynamics 
along tilted periodic potentials was recovered by setting 
f b (X,b t ) = -V Xz [(V /2)(l-cos(2irx l /\))-b t x l }, and we 
used the same parameters as in the single-parameter re- 
sampling case. We found that using 10 6 trajectories the 
discrepancies with respect to the data in Table II were 
about 16 % at largest. This is due to the fact that when 
simulations are run on a smooth surface, all trajecto- 
ries of the dimer are weighted equally. However, in the 
actual periodic potential the main contribution to the 
diffusion coefficient comes from paths crossing the sad- 
dle point from one minimum to another. Thus, in the 
resampling procedure most of the paths are not relevant 
for determining the value of D, and thus the errors re- 
main relatively large even with 10 6 paths for the present 
case. 



SUMMARY AND CONCLUSIONS 

In this work, we have employed the recently proposed 
SPHD scheme to study the diffusive motion of Brown- 
ian particles in periodic potentials in ID. Unlike the HD 
schemes proposed so far, the SPHD scheme allows an ex- 
act correction of the biased dynamics based on reweigh- 
ing of all the transition paths. There's also no restric- 
tion for the type of bias potential used, and thus the 
SPHD method can be used to extrapolate results to mul- 
tiple values of the bias force from a single simulation run, 
as shown by the two cases we have studied here: time- 
dependent forcing for a single Brownian particle and con- 
stant forcing for a Brownian dimer. For a Brownian par- 
ticle in an external ac bias, our results are in agreement 



with previous studies and show that there's no stochastic 
resonance in this system. For the second case of a dimer 
moving in a periodic potential with a static bias we find 
a strong suppression of D when the dimer bond length is 
an integer multiple of the lattice constant of the poten- 
tial. This suppression is weakened by an applied dc bias. 
Our work demonstrates how external forcing can be used 
to control particle mobilities in periodic potentials. 
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